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EXECUTIVE  SUMMARY 


Our  3-dimensional,  time-dependent  model  of  the  global  ionosphere  was  streamlined  in 
an  effort  to  develop  a  computationally  fast,  user  Mendly,  reliable,  Ionospheric  Forecast 
Model  (IFM)  for  Air  Weather  Services.  The  model  yields  predictions  for  the  NO,  O2+, 
N2+,  and  0+  densities  and  the  ion  and  electron  temperatures  over  the  globe  at  E  and  F 

region  altitudes.  The  model  also  contains  a  simple  algorithm  for  predicting  H+  densities  in 
the  F  region.  The  inputs  needed  by  the  IFM  are  global  distributions  of  the  neutral 
densities,  temperatures  and  winds,  the  auroral  oval  precipitation,  the  magnetospheric  and 
dynamo  electric  fields,  and  the  topside  electron  heat  flux.  Because  of  the  model's  modular 
construction,  the  IFM  can  readily  accept  different  global  input  patterns  and,  hence,  has  the 
capability  of  being  driven  by  retd  -  time  inputs  from  Air  Force  satellites  or  ground  -  based 
sites.  In  the  current  version  of  the  model,  the  input  patterns  have  been  selected  and  the 
IFM  is  therefore  self  -  contained  and  can  be  driven  by  simple  geophysical  indices.  With  a 
simple  specification,  the  IFM  ou^uts  a  variety  of  parameters,  including  global  snapshots  of 
NinF2  and  hmF2,  density  and  temperature  contours  at  fixed  altitudes,  altitude  profiles  at 
given  locations,  as  well  as  other  parameters. 
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1.  INTRODUCTION 


Over  a  20  year  period,  we  have  been  actively  involved  in  developing  numerical  models 
of  the  terrestrial  ionosphere.  The  most  comprehensive  model  that  we  constructed  was  a 
time  -  dependent,  three  -  dimensional,  multi  -  ion  model  of  the  global  ionosphere  at  E  and 
F  region  altitudes  (cf.  Schunk,  1988a;  Sojka,  1989).  This  model  yields  density 

distributions  for  six  ion  species  (NO^■,  O2+,  N2+,  0+,  N+,  He+)  and  Ae  electron  and  ion 
temperature  distributions  via  a  numerical  solution  of  the  coupled  continuity,  momentum, 
and  energy  equations.  The  model  takes  account  of  field-aligned  diffusion,  cross-field 
electrodynamic  drifts  both  in  the  equatorial  region  and  at  high  latitudes,  interhemispheric 
flow,  thermospheric  winds,  polar  wind  escape,  energy-dependent  cheimcal  reactions, 
neutral  composition  changes,  ion  production  due  to  solar  EUV  radiation,  resonantly 
scattered  radiation  and  auroral  precipitation,  electric  field  induced  temperature  aiusotropies, 
thermal  conduction,  diffusion-thermal  heat  flow,  and  local  heating  and  cooling  processes. 
The  model  also  takes  account  of  the  offset  between  the  geomagnetic  and  geo^aphic  poles 
in  both  the  northern  and  southern  hemispheres  and  the  bending  of  field  lines  near  the 
magnetic  equator. 

With  the  funding  from  this  Air  Force  contract,  we  streamlined  the  above  global 
ionospheric  model  in  order  to  construct  a  conaputationally  efficient,  user  friendly. 
Ionospheric  Forecast  Model  for  Air  Weather  Services.  The  development  effort  involved 
four  major  tasks.  The  first  task  was  primarily  concerned  with  the  development  of  the 
ionospheric  code  for  a  single  convecting  flux  tube.  This  part  of  the  code  had  to  redone 
so  that  it  could  be  made  computationaiUy  fast  and  efficient.  This  involved  deriving  the 
ambipolar  diffusion  equations  for  two  ions,  modifs^g  the  ion  chemistry,  and  simplifying 
the  ion  and  electron  energy  equations.  The  numerical  inversion  procedure  also  had  to  be 
redone. 

The  second  task  involved  the  development  of  the  codes  that  describe  the  inputs  to  the 
ionospheric  model  as  well  as  the  effort  involved  in  linking  the  single  flux  tube  code  to  the 
coordinate  systems  and  global  inputs.  Decisions  had  to  be  made  about  what  coordinate 
systems  to  adopt  and  how  many  flux  tubes  to  follow.  The  inputs  to  the  ionospheric 
forecast  model  that  needed  to  be  considered  were  a  neutral  wind  pattern,  a  neutral 
atmosphere,  magnetospheric  and  equatorial  electric  field  distributions,  auroral  particle 
precipitation  patterns,  the  topside  electron  heat  flux,  and  the  mid-latitude  protonospheric 
exchange  flux.  Although  many  of  these  inputs  were  readily  available,  they  had  to  be 
modified  so  as  to  make  them  computationally  efficient 

The  third  task  involved  an  extensive  series  of  validation  runs.  The  Ionospheric 
Forecast  Model  was  tested  against  most  of  the  54  diumdly  reproducible  runs  already 
conducted  with  our  full  ionospheric  model  for  the  High  Latitude  Ionospheric  Specification 
Model  (HLISM).  The  IFM  was  also  vaUdated  as  pm  of  the  PRIMO  ^d  VIM  programs. 
The  last  task,  involved  the  development  of  graphic  displays  and  user  friendly  I/O  interface 
software. 

In  the  following  sections  of  this  report,  we  provide  a  complete  description  of  the 
current  version  of  the  IFM.  We  also  briefly  indicate  the  operating  mode  of  the  model,  but 
a  user  manual  giving  detailed  operating  instructions  will  be  provided  separately. 
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2.  STREAMLINING 


As  indicated  earlier,  a  significant  effort  was  devoted  to  streamlining  our  global 
ionospheric  models,  and  the  details  are  provided  in  the  subsections  that  follow. 


2..1.  Piffusigg  Equations 

The  mathematical  formulation  of  our  'comprehensive'  ionospheric  model  is  based  on 
the  continuity,  momentum,  and  energy  equations  for  four  major  ions  (NO+,  O2+,  N2+, 

0+)  and  two  minor  ions  (He+,  N+).  The  ion  diffusion  equations  are  obtained  from  the 
momentum  equations  by  assuming  that  the  plasma  flow  is  ambipolar  and  subsonic,  which 
are  valid  assumptions  at  E  and  F  region  altitudes.  In  the  original  formulation,  the  four 
coupled  diffusion  equations  are  solved  simultaneously  over  the  complete  altitude  range 
from  100  to  800  km. 

In  the  streamlined  model,  the  E  and  F  regions  of  the  ionosphere  are  formulated 
separately.  In  the  E  region,  we  take  account  of  the  fact  that  dr^usion  is  negligible  and 
obtain  the  ion  densities  by  assuming  that  chemical  equilibrium  prevails, 

ni  =  Pi/Li  (1) 

where  Pi  is  the  ion  production  rate,  Li  the  ion  loss  frequency,  and  ni  the  ion  density.  Here, 
the  two  minor  ions  He+  and  N+  are  neglected,  but  the  remaining  four  ions  (NO+,  02‘’', 
N2+,  0+)  are  treated  as  major  ions  (i.e.,  equally  important). 

In  the  F  region,  the  streamlining  effort  centered  around  the  assumption  that  only  NO+ 

and  0+  are  major  ions,  which  has  been  verified  over  the  years  via  numerous  runs  of  our 
comprehensive  ionospheric  model.  Therefore,  the  original  triple  -  ion  diffusion 
formulation  can  be  replaced  with  the  following  reduced  two  -  ion  formulation; 


|^  +  |(n,„,)  =  P,-L,n, 

Bn.  d  ^  X  „  , 

(n2U2 )  —  P2  ~  L2n2 

Uj  =  Vj^  -  D,  sin^  ] 

Uj  =  -02^*^^  I 

He  =  ni  +  n2 
ngUe  =  niui  +  n2U2 


L^  +  iM+±±(T  +T  .  I  (T./T,)3n. 
Ln,  az  kT,  T.az' ■  IT 


1  dn 


n^  dz 


kT 


T2  dz 


dz 


(2) 

(3) 

(4) 

(5) 
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where  the  subscripts  1,2,  and  e  refer  to  NO+,  O  and  electrons,  respectively,  and  where  ui 
is  the  vertical  ion  drift  velocity,  mi  is  the  ion  mass,  Ti  is  the  ion  temperature,  Te  is  the 

electron  temperature,  Di  =  kTi/(mi\)i)  is  the  ion  diffusion  coefficient,  Viz  is  the  induced 

vertical  ion  drift  due  to  neutral  winds  and  electric  fields,  \)i  =  ZnVin  is  the  total  ion  collision 
frequency  with  neutrals,  k  is  Boltzmann's  constant,  g  is  the  gravitational  acceleration,  z  is 
the  vertical  coordinate,  t  is  time,  and  I  is  the  inclination  angle  of  the  geomagnetic  field. 

In  the  streamlined  model,  the  minor  ions  N+  and  He+  are  ignored  and  the  two 
remaining  molecular  ions  (O2+,  N2+)  are  taken  to  be  minor  ions  in  chemical  equilibrium. 
Therefore,  their  densities  are  obtain^  from  equation  (1)  at  all  altitudes  throughout  the  E 
and  F  regions. 

In  addition  to  the  reduction  of  the  diffusion  formulation  from  a  triple-ion  to  a  two-ion 
scheme,  several  other  simplifications  were  made.  The  small  effect  of  temperature 
anisotropies  was  neglected,  and  only  ion  collisions  with  the  dominant  neutrals  (N2,  O) 
were  included  in  the  total  ion  collision  frequency. 

Note  that  the  above  diffusion  formulation  is  valid  only  at  middle  and  high  latitudes, 
where  the  geomagnetic  field  lines  are  straight  but  inclined  to  the  vertical.  At  equatorial 
latitudes  a  different  formulation  is  employed,  which  will  be  discussed  later.  Also  note  that 
the  separation  of  the  E  and  F  region  formulations  ^ows  for  different  numerical  step  sizes 
in  the  two  domains,  which  has  a  distinct  computational  advantage.  This  will  be  discussed 
in  more  detail  later. 


2.2.  Ion  Chemistry 

The  ion  chemistry  in  our  comprehensive  ionospheric  model  is  fairly  involved.  The 
chemical  scheme  includes  numerous  energy  -  dependent  ion  -neutral  chemical  reactions, 
electron  -  ion  recombination,  photoionization,  ionization  due  to  auroral  electron 
precipitation  and  ionization  resulting  from  resonantly  scattered  solar  radiation.  In  the  IFM, 
these  processes  have  been  simplified  and/or  improv^. 

Because  some  minor  ions  were  neglected  in  the  construction  of  the  IFM,  the  ion  - 
neutral  chemical  reaction  scheme  could  be  shortened.  Also,  over  the  years  we  found  that 
some  of  the  chemical  reactions  included  in  the  comprehensive  ionospheric  model  were,  in 
fact,  negligible.  Therefore,  the  IFM  could  be  constructed  with  a  reduced  ion  -  neutral 
chemical  scheme,  which  is  given  in  Table  1. 

With  regard  to  photoionization,  rigorous  photoproduction  calculations  are  performed 
over  the  entire  altitude  range  from  100  to  800  km  in  the  comprehensive  ionospheric  model. 

Also,  these  calculations  are  performed  even  for  large  solar  zenith  angles  (x  >  90o),  which 
requires  complex  expressions  for  grazing  incidence  functions.  However,  experience  has 
shown  that  important  simplifications  can  be  made  without  loss  of  accuracy.  First,  above 
300  km  photoionization  frequencies  can  be  introduced.  The  ion  photoionization  frequenaes 
are  calculated  rigorously  at  300  km  and  then  these  values  are  used  to  calculate  ion 
production  rates  above  this  altitude.  Also,  ion  production  rates  are  only  calculated  for  solar 

zenith  angles  less  than  90°,  which  eliminates  the  need  for  the  complex  grazing  incidence 
functions. 
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The  calculation  of  the  photoionization  rates  requires  a  knowledge  of  the  solar  spectrum 
and  the  absorption  and  ionization  cross  sections  for  the  neutral  gases  in  the  earth's  upper 
atmosphere.  The  values  adopted  for  the  IFM  were  obtained  from  Torr  et  al  (1985),  but  the 
37-  wavelength  intervals  were  adjusted  to  11  intervals.  However,  checks  were  made  to 
verify  that  the  production  rates  obtained  from  the  1 1-  interval  scheme  were  within  a  few 
percent  of  those  obtained  from  the  37  intervals.  Also,  in  the  BFM,  the  solar  fluxes  are 
scaled  with  Fio.7  in  order  to  account  of  solar  cycle  effects. 

In  our  comprehensive  ionospheric  model,  the  ion  production  rates  due  to  auroral 
electron  precipitation  are  obtained  by  scaling  the  production  rate  profiles  obtained  from 
Knudsen  et  al  (1977)  with  the  auroral  electron  energy  flux.  In  the  IFM,  our  calculation  of 
auroral  production  rates  has  been  improved  significantly.  Our  new  production  rate  profiles 
are  based  on  the  auroral  electron  deposition  code  developed  by  Strictland  et  al  (1976). 
R.E.  DanieU  of  CPI  ran  this  code  for  numerous  cases  and  then  gave  us  the  resulting  sets  of 
production  rate  proves.  Specifically,  he  considered  27  geophysical  situations,  covering  3 
levels  of  solar  activity,  3  seasons,  and  3  levels  of  geomagnetic  activity.  For  each  of  the  27 
geophysical  situations  he  considered  5  auroral  cases,  covering  different  characteristic 
energies.  For  the  135  cases,  production  rate  profiles  were  obtained  for  N2+,  O2+  and  0+, 
and  we  incorporated  these  profiles  in  a  subroutine  for  efficient  access  when  plasma  flux 
tubes  convect  through  the  auroral  oval  and  encounter  changing  auroral  conditions. 

In  the  IFM,  the  nocturnal  production  rates  due  to  resonantly  scattered  solar  radiation  are 
identical  to  those  contained  in  our  comprehensive  ionospheric  model.  These  rates  are  those 
presented  by  Chen  and  Harris  (1971),  which  do  not  vary  during  the  lught. 


2.3.  Ion  Energy  Equation 

The  ion  temperature  calculation  in  our  comprehensive  ionospheric  model  is  based  on  a 
numerical  solution  of  a  nonlinear,  second  order,  partial  differential  equation.  Such  an 
equation  results  from  the  inclusion  of  thermal  conduction  and  diffusion  -  thermal  transport 
processes.  However,  these  transport  processes  only  start  to  have  an  effect  at  altitudes 
above  600-700  km  and  only  under  certain  circumstances.  Since  an  appreciable  CPU 
savings  can  be  obtained  by  neglecting  these  processes,  they  are  not  included  in  the  IFM. 
In  this  case,  the  ion  energy  equation  reduces  to  a  balance  between  local  heating  (Qi)  and 
cooling  (Ci)  processes. 


Qi  =  Ci  (8) 

As  with  our  comprehensive  ionospheric  model,  the  calculation  of  the  ion  heating  and 
cooling  rates  in  the  IFM  includes  thermal  coupling  to  the  electrons  and  neutrals  as  well  as 
electric  field  heating.  However,  some  simplifications  were  possible.  Only  coupling  to  the 
dominant  neutrals  N2  and  O  was  included  and  the  small  effect  associated  with  ion 
temperature  anisotropies  was  neglected. 

2.4.  Electron  Energy  Equation 

Thermal  conduction  is  a  very  important  process  for  the  electron  gas  at  F  -  region 
altitudes  and  cannot  be  neglected.  Therefore,  in  contrast  to  the  Ti  situation,  a  nonlinear, 
second  order,  parabolic,  partial  differential  equation  must  be  solved  in  order  to  obtain 
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electron  ten^ratures.  This  equation  can  be  expressed  in  the  form 


+  Q.-Ce 


(9) 


where  Ke  is  the  electron  thermal  conductivity,  Qe  is  the  electron  heating  rate,  and  Ce  is  the 
electron  cooling  rate.  In  addition  to  thermal  conduction,  our  comprehensive  ionospheric 
model  included  electron  heating  via  both  auroral  electrons  and  photoelectrons,  thermal 
coupling  to  the  ions,  cooling  via  elastic  collision  with  five  neutral  species  (N2,  O2,  O,  He, 
H),  and  cooling  via  several  inelastic  collisional  processes  (N2  &  O2  rotational  excitation, 
N2  &  O2  vibrational  excitation,  and  excitation  of  the  fine  structure  levels  of  atomic 
oxygen). 

In  the  IFM,  only  the  dominant  cooling  processes  are  retained.  These  are  elastic 
collisions  with  N2  and  O,  and  excitation  of  the  fine  structure  levels  of  atomic  oxygen.  The 
neglect  of  the  other  unimportant  cooling  terms  results  in  a  significant  CPU  savings  because 
some  of  these  terms  are  long  and  complex.  The  expressions  for  the  processes  included  are 
given  in  Schunk  (1988a)  and  are  not  rq)eated  here. 


2.5.  Equatorial  Model 

The  mathematical  formulation  described  above  is  restricted  to  middle  and  high  latitudes 
because  it  is  based  on  the  assumptions  that  vertical  ^adients  dominate  and  that  the  B-field 
lines  are  straight  but  inclined  to  the  vertical.  In  this  case,  the  plasma  transport  equations 
can  be  solved  as  a  function  of  altitude  for  convecting  plasma  flux  tubes.  At  equatorial 
latitudes,  on  the  other  hand,  the  geomagnetic  field  varies  from  being  inclined  to  the  vertical 
at  higher  latitudes  to  being  horizontal  at  the  equator.  Therefore,  the  transport  equations 
must  be  solved  along  the  field  lines  taking  into  account  both  horizontal  and  vertical 
gradients.  In  our  comprehensive  model  of  the  global  ionosphere,  we  adopted  the  Sterling 
et  al  (1969)  equatorial  ionospheric  model  and  modified  the  diffusion  coefficients  and  ion 
chemistry  so  that  they  were  consistent  with  those  used  in  our  mid-high  latitude  ionosphere 
model  (cf.  Sojka  and  Schunk,  1985). 

In  the  IFM,  we  replaced  the  Sterling  et  al  (1969)  model  with  the  Anderson  (1973, 
1981)  equatorial  ionospheric  model.  With  this  model,  the  continuity  and  momentum 

equations  for  0+  and  electrons  are  solved  along  a  magnetic  field  line  that  extends  fi:om  the 
northern  to  the  southern  hemisphere.  The  model  takes  into  account  ion-neutral  chemical 
reactions,  photoionization,  diffusion,  neutral  winds,  and  electrodynamic  drifts.  In  the 
most  recent  version,  the  model  automatically  considers  a  series  of  plasma  flux  tubes  with 
different  equatorial  crossing  altitudes  so  that  altitude  profiles  can  be  obtained  at  any 
location.  However,  the  series  of  flux  tubes  is  at  only  one  longitude.  Typically,  the 
equations  are  solved  as  a  function  of  time  along  the  series  of  flux  tubes  as  they  corotate  and 
convect  in  response  to  dynamo  electric  fields.  The  lower  boundary  conditions  are 
generally  taken  to  be  at  90  km  in  both  hemispheres. 

When  the  Anderson  (1973)  equatorial  ionospheric  model  was  incorporated  into  the 
IFM,  several  modifications  were  necessary.  First,  our  E-region  model  is  global  and  covers 
the  altitude  range  fi'om  90-160  km,  so  the  lower  boundary  in  the  Anderson  model  was 
moved  up  to  160  km.  Next,  the  photochemistry  and  diffusion  coefficients  were  modified  to 
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be  consistent  with  those  used  in  the  mid  -  high  latitude  part  of  the  IFM.  Also,  the  model 
was  extended  from  one  longitude  to  aU  longitudes  (global)  and  the  input  routines  were 
modified  to  accept  time-dependent,  global  inputs. 


Table  1  Ion  Chemistry  and  Reaction  Rates 


Reaction 

Rate  Coefficients,  cm3  s-i 

0+  +  N2  -4  N0+  +  0 

k\ 

O'**  +  O2  O2***  0 

ki 

O2+  +  N2  -4  N0+  +  NO 

5  X  l(hi6 

^2^  +  €  — >  0  +  0 

1.6  X  10-7  (300/1^)0.55 

N2+  +  0  ->  N0+  +  N 

1.4  X  10-10  (300/7^0.44;  j<  1500*  K 
5.2  X  10-11  (7/300)0.2;  T>  1500*  K 

N2+  +  0  ->  0+  +  N2 

1  X  10-11  (300/7)0.23;  T<  1500*  K 

3.6  X  10-12  (7/300)0.41;  7>  1500*  K 

N2+  +  O2  -4  O2+  +  N2 

5  X  10-11  (300/7) 

N2+  +  c  ^  N  +  N 

1.8  X  10-7  (300/7^)0.39 

N0+  +  e  -4  N  +  0 

4.2  X  10-7  (300/7^)0.85 
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3.  IFM  INPUTS 


The  IFM,  like  all  ionospheric  models,  requires  both  magnetospheric  and  thermospheric 
inputs  in  order  to  account  for  ionospheric  coupling  to  these  domains.  The  inputs  need  to  be 
time-dependent  and  global.  At  present,  various  empirical  (i.e.,  statistical)  models  are 
available  for  almost  ^1  of  the  required  inputs,  and  these  input  models  can  be  'adjusted'  to 
accommodate  real-time  measurements  made  by  Air  Force  satellites  or  at  Air  Force  ground 
sites.  Since  it  was  anticipated  that  significant  in^rovements  in  the  inputs  would  be  made  in 
the  coming  years,  the  IFM  was  constructed  in  a  modular  form  so  that  different  input 
models  could  be  'plugged  in'  with  no  impact  on  the  rest  of  the  IFM.  However,  for  the 
current  version  of  the  IFM,  a  specific  set  of  input  models  was  adopted.  This  first  version 
of  the  IFM  is  therefore  self-contained  and  dnven  by  simple  geophysical  indices.  The 
specific  input  models  adopted  are  briefly  described  in  the  following  subsections. 


3.1.  Neutral  Atmosphere 

The  IFM  requires  a  global  distribution  of  the  neutral  densities  (N2,  O2,  O)  and 
temperature  at  altitudes  between  100-800  km.  These  input  parameters  are  obtained  from 
the  MSIS-90  atmospheric  model  (Hedin,  1991).  This  empirical  model  yields  the 
atmospheric  parameters  for  different  solar  cycle,  seasonal,  and  geomagnetic  activity  levels. 
It  also  contains  diurnal  and  longitudinal  variations. 


3.2.  Neutral  Wind 

The  neutral  wind  model  adopted  as  an  input  to  the  IFM  is  the  MSIS-Wind  model 
constructed  by  Hedin  et  al  (1991).  This  empirical  model  provides  global  distributions  of 
the  zonal  and  meridional  wind  components.  As  with  the  atmospheric  model,  the  wind 
model  takes  account  of  solar  cycle,  seasonal,  geomagnetic  activity,  diurnal,  and 
longitudinal  variations. 


3.3.  Magnetospheric  Electric  Field 

At  high  latitudes,  the  plasma  convection  pattern  induced  by  magnetospheric  electric 
fields  is  crucial  in  determining  the  ionospheric  densities  and  temperatures.  Because  of  its 
importance,  several  empirical  models  of  the  convection  electric  field  have  been  developed 
over  the  years.  The  m^els  have  been  based  on  both  radar  data  and  satellite  electric  field 
data  (Volland,  1978;  Foster,  1983;  Heelis  et  al,  1982;  Heelis,1984;  Heppner  and 
Maynard,  1987).  For  the  current  version  of  IFM,  we  adopted  the  Heppner  and  Maynard 
1987)  electric  field  model.  This  model  provides  magnetospheric  electrostatic  potential 
distributions  over  the  polar  regions  as  a  function  of  Kp  for  southward  IMF  (Bz<0)  and 
both  negative  and  positive  By  values. 


3.4.  Auroral  Precipitation 

Particle  precipitation  in  the  auroral  oval  acts  as  an  ionization  source,  a  source  of  bulk 
heating  for  the  electron  gas,  and  a  thermal  conduction  source  through  our  upper  boundary. 
As  with  the  electric  field,  several  empirical  models  are  available  that  describe  the 
precipitating  electron  energy  flux  and  characteristic  energy  (Feldstein  and  Starkov,  1967; 
Spiro  et  al,  1982;  Wallis  and  Budzinski,  1981;  Hardy  et  al,1985).  For  the  IFM,  we 
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adopted  the  model  by  Hardy  et  al  (1985),  which  is  the  most  comprehensive  model  of 
electron  precipitation  developed  to  date.  In  the  IFM,  as  the  plasma  flux  tubes  move 
through  the  auroral  oval,  the  characteristic  energy  and  energy  flux  values  obtained  from  the 
Hardy  et  al  (1985)  model  are  used  to  calculate  particle  production  rates,  the  bulk  electron 
heating  rate,  and  the  thermal  conduction  flux  ^ough  our  upper  boundary  at  the  current 
location  of  the  flux  tube. 


3.5.  Equatorial  Electric  Field 

Dynamo  electric  fields,  generated  via  the  thermospheric  wind,  are  extremely  important 
for  determining  equatorial  electron  densities.  Until  recently,  the  best  empirical  model  of 
equatorial  electric  fields  in  the  one  developed  by  Richmond  et  al  (1980),  which  is  the  one 
used  in  our  comprehensive  ionospheric  model.  Recently,  however,  B.  Fejer  and 
colleagues  constmcted  a  sophisticated  empirical  model  of  equatorial  electric  fields  including 
longitudinal  variations,  and  the  ionospheric  modeUing  results  obtained  with  these  electric 
fields  have  been  impressive  (Preble  et  al,  1994).  Therefore,  the  Fejer  equatorial  electric 
field  model  was  adopted  for  the  IFM. 


3.6.  Topside  Heat  Flow 

Escaping  photoelectrons  and  auroral  electrons  lose  energy  to  the  thermal  electrons  at 
high  altitudes  as  they  traverse  the  ionospheric  plasma.  The  transferred  energy  is  then 
conducted  down  to  the  lower  ionosphere  and  enters  as  a  topside  (800  km)  boundary 
condition  in  the  Ionospheric  Forecast  Model.  The  procedure  for  calculating  the  boundary 
heat  flow  values  is  identical  to  that  used  in  our  comprehensive  ionospheric  model  and  is  not 
repeated  here  (cf.  Schunk  et  al,  1986). 
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4.  COORDINATE  SYSTEMS  AND  TRAJECTORIES 


The  Ionospheric  Forecast  Model  is  a  Lagrange  -  Euler  hybrid  model  in  that  the 
continuity,  momentum,  and  energy  equations  are  solved  on  a  fixed  spatial  grid  either  as  a 
function  of  altitude  (mid  and  high  latitudes)  or  along  geomagnetic  field  lines  (low  latitudes) 
for  £  X  B  convecting  flux  tubes  of  plasma  (cf.  Figure  1).  As  a  consequence,  coordinate 
systems  and  trajectory  selection  are  in^ortant  aspects  of  the  IFM. 


4.1.  Coordinate  Systems 

There  are  two  primary  coordinate  systems  associated  with  the  IFM,  namely  geographic 
and  geomagnetic.  Our  comprehensive  ionospheric  model  has  always  been  bas«I  on  a  tilted 
offset  dipole  magnetic  system  because  relatively  simple  transformations  exist  to  convert 
firom  this  frame  to  geographic.  Also,  we  have  previously  assumed  that  the  dipole  magnetic 
coordinate  system  is  equivalent  to  the  magnetic  invarient  system.  In  the  past,  this 
assumption  was  justified  because  most  of  our  studies  were  parametric  in  nature,  and  with 
latitude  bins  of  3  degrees  and  local  time  bins  of  1  hour,  the  errors  associated  with  this 
assumption  were  not  significant.  However,  in  future  applications  with  real-time  data 
(images,  in  situ  satellite  measurements  and  ground-based  observations),  this  assumption  is 
inappropriate.  Hence,  in  the  IFM,  the  more  rigorous  IGRF  magnetic  field  model  is 
adopted,  and  it  is  represented  by  bicubic  spline  interpolation  tables  for  the  magnitude, 
inclination,  and  declination.  The  CGMC  transformation  used  at  Phillips 
Laboratory/Boston  is  adopted  for  transforming  between  the  geographic  and  magnetic 
frames. 


4.2.  Trajectory  Selection 

The  magnetospheric  electric  field  models  provide  contours  of  electrostatic  potential, 
which  coincide  with  the  plasma  streamlines  at  F-region  altitudes  (cf.  Figure  1).  However, 
decisions  have  to  be  made  with  regard  to  how  many  trajectories  to  select,  how  many 
plasma  flux  tubes  to  follow,  and  where  the  flux  tubes  should  be  placed  initially.  Also, 
account  must  be  taken  of  the  fact  that  as  the  plasma  flux  tubes  convect  in  response  to  the 
imposed  magnetospheric  electric  field  pattern,  their  spatial  distribution  will  become  non- 
un&brm  because  the  various  trajectories  have  different  lengths  and  because  the  magnitude 
of  the  electric  field  varies  around  each  trajectory  differently.  Therefore,  a  multi-step 
process  is  needed  to  set  up  and  follow  convecting  plasma  flux  tubes. 

In  the  IFM,  the  procedure  used  is  shown  in  Figure  2.  First,  a  latitude-longitude  grid  is 
initialized  (INTTGRID).  Then,  the  Kp  and  other  proxy  indices  are  acquired  as  a  time  series 
from  AWS  data  computers  (INITACTV).  These  data  are  then  used  by  INITTRAJ  to 
initialize  the  plasma  flux  tubes  diat  are  to  be  followed.  Flux  tubes  are  placed  at  the  center 
of  grid  cells.  FLOW  generates  the  trajectories  that  the  plasma  flux  tubes  are  to  follow. 
n^OW  calls  FLOWVl  to  obtain  the  corotational  electric  field  and  FLOWV3  to  obtain  the 
Heppner-Maynard  electric  field  pattern.  As  the  flux  tubes  convect  in  response  to  the 
imposed  electric  fields,  their  spatial  distribution  can  change  markedly.  Hence,  eve^  15 
minutes  of  model  time,  it  is  necessary  to  check  the  spatial  distribution  of  flux  tubes  in  the 
ionospheric  grid  to  m^e  sure  there  is  adequate  coverage  everywhere.  This  is  done  in 
CELLCHECK.  Flux  tube-  are  deleted  if  too  many  appear  in  a  spatial  cell  in  order  to  save 
CPU  time.  Flux  tube  are  added  to  a  spatial  cell  tf  there  is  an  insufficient  number  of  flux 
tube  in  that  cell.  This  is  done  in  steps.  FILLCELL  creates  a  new  trajectory  that  passes 
through  the  center  of  the  empty  cell  and  obtains  the  coordinates  (CELLCOORDS).  The 
trajectory  is  followed  backwards  in  time  to  the  previous  model  time  interval  and 
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interpolation  is  used  to  obtain  the  flux  tube  density  profile  at  that  previous  time  interval 
(FLOWBACK).  This  flux  tube  reverses  direction  and  then  follows  the  trajectory  back  to 
the  present  time  (REVERSE).  Of  course,  it  ends  up  at  the  center  of  the  empty  cell.  The 
routine  FILLHIST  creates  a  record  of  the  distribution  of  plasma  flux  tubes  after  each  IFM 
time  interval  for  diagnostic  purposes.  Finally,  before  following  the  trajectories  to  the  next 
model  time  interval,  the  model  grid  is  advanced  in  MLT  by  the  time  interval 
(ANOTHERFLOW).  This  assures  that  corotational  flux  tubes  stay  at  the  midpoints  of  their 
cells  and,  thereby,  reduces  the  number  of  flux  tubes  needed. 

An  example  of  how  the  flux  tube  selection  procedure  works  is  shown  in  Figures  3  -  5. 
Initially,  a  spatial  grid  with  cells  that  are  5o  in  latitude  and  1  hour  in  MLT  is  created  for  the 

region  poleward  of  40°  latitude.  One  plasma  flux  tube  is  placed  at  the  center  of  each  of  the 
120  cells  (Figure  3,  bottom  panel).  Subsequently,  the  flux  tubes  convect  in  response  to  the 
corotational  and  magnetospheric  electric  fields.  In  this  test  case,  a  Volland  (1978) 
symmetric  two  -  ceU  convection  pattern  is  used.  As  time  evolves,  the  distribution  of  flux 
tubes  becomes  nonuniform.  At  15-minute  time  intervals,  the  distribution  is  checked  and 
cells  with  10  or  more  flux  tubes  have  the  number  reduced  to  2-5,  while  cells  with  an 
insufficient  number  are  augmented  with  new  flux  tubes.  The  resulting  flux  tube 
distribution  after  1  day  is  shown  in  Figure  3  (middle  panel)  and  after  3  days  in  Figure  3 
(top  panel).  It  is  clear  that  there  is  always  adequate  coverage  of  the  high  latitude  region  for 
this  Kp  =  1  Volland  convection  pattern.  Similar  results  are  shown  in  Figures  4  and  5  for 
Kp  =  4  and  7,  respectively.  Again,  our  flux  tube  selection  procedure  provides  adequate 
coverage  over  the  entire  polar  region.  Note  that  the  maximum  number  of  flux  tubes  in  400 
for  Kp  =  7. 
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Figure  1.  A  schematic  view  of  the  F-region  ionosphere  bounded  below  (100km)  and 
above  (800km)  extending  from  the  north  pole  to  the  magnetic  equator.  T^e  imposed 
magnetospheric  electric  fields  and  auroral  oval  are  shown  as  contoured  lines  and  shaded 
regions,  respectively.  A  midnight  cross  section  reveals  the  direction  of  motion  of  a  high 
latitude  (Vi),  mid-latitude  (V2),  and  equatorial  (Vs)  plasma  flux  tube.  From  Sojka  [1989]. 


Figure  2.  Block  diagram  showing  how  many  trajectories  to  select,  how  many  plasma 
flux  tubes  to  follow,  and  where  the  flux  tubes  should  be  placed  initially. 
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Figure  3.  Flux  tube  locations  at  the  start  of  the  simulation  (bottom  panel),  after  1  day 
(middle  panel),  and  after  3  days  (top  panel)  for  Kp  =  1. 
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Figure  4.  Flux  tube  locations  at  the  start  of  the  simulation  (bottom  panel),  after  1  day 
(middle  panel),  and  after  3  days  (top  panel)  for  Kp  =  4. 
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Figure  5.  Flux  tube  locations  at  the  start  of  the  simulation  (bottom  panel),  after  1  day 
(middle  panel),  and  after  3  days  (top  panel)  for  Kp  =  7. 
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5.  NUMERICAL  SOLUTION 


For  the  mid-lugh  latitude  domain,  the  diffusion  and  thermal  conduction  equations  are 
solved  as  a  function  of  altitude  for  convecting  plasma  flux  tubes.  In  our  comprehensive 
ionospheric  model,  a  fixed  spatial  step  is  used  from  the  bottom  to  the  top  boundary  (100- 

800  km).  Also,  three  ions  (NO,  O2+,  O)  are  assumed  to  be  major  ions,  which  means 
the  three  coupled  diffusion  equations  must  be  solved  simultaneously.  However, 

experience  has  shown  that  NO  and  O  are  the  dominant  ions  in  the  altitude  region  where 
diffusion  is  important,  and  a  2-ion  code  is  computationally  more  efficient  than  a  3-ion 
code.  Hence,  the  IFM  is  based  on  a  two  "major  ion"  formulation  in  the  F-region.  Also,  a 
small  spatial  step  is  generally  needed  in  the  E-region  to  resolve  the  peak,  but  a  larger  spatial 
step  can  be  used  in  the  F-region  because  the  peak  is  much  broader.  Again,  for 
computational  efficiency,  the  IFM  formulation  is  divided  into  separate  E  and  F  regions. 

The  IFM  solution  procedure  is  as  follows.  An  altitude  is  selected  which  defines  the 
separation  between  the  E  and  F  regions  ( «130  km).  In  the  E-region  a  4  km  spatial  step  is 
used,  while  in  the  F-region  the  spatial  step  is  larger  (10-20  km),  as  shown  in  Figure  6, 
The  E  and  F  region  density  ^uations  are  solved  as  a  function  of  time  with  a  10  -  100 
second  time  step  for  convecting  flux  tubes.  Every  15  minutes,  the  required  inputs,  the 
associated  parameters,  and  the  plasma  temperatures  are  calculated.  Tlie  sequence  of  these 
calculations  is  shown  in  Figme  7.  Subroutine  ONE  first  calls  routines  that  calculate 
altitude  profiles  for  the  neutral  densities  and  temperature  (ATMOS),  the  electron 
temperature  (TE),  the  ion  temperature  (TI),  diffusion  coefficients  (DIFF),  chemical 
reactions  (CKCEM),  photoionization  (PHOTO  and  SOLAR),  production  due  to  resonantly 
scattered  solar  radiation  (RESONAh^  and  auroral  precipitation  (AURORA),  the  bottom 
boundary  conditions  (BC),  and  the  wind-induced  vertical  ion  drift  (WIND).  With  altitude 
profiles  of  these  parameters,  ONE  then  calculates  the  coefficients  that  appear  in  the  partial 
differential  equations  that  describe  the  spatial  and  temporal  variation  of  the  ion  densities. 


5.1.  E-Region 

In  the  above  solution  procedure,  the  E-region  densities  are  obtained  by  assuming 
chemical  equilibrium  (equation  1),  but  all  four  ions  are  allowed  to  have  comparable 
densities.  The  resulting  four  coupled  nonlinear  equations  are  solved  by  first  expanding  the 
coupling  and  nonlinear  terms  in  a  Taylor  series  in  time  and  then  by  iterating  to  a  solution 
from  an  initial  guess.  The  iteration  can  be  done  separately  at  each  altitude  step. 


5.2.  F-Region 

The  F-region  densities  are  obtained  from  diffusion  equations  for  the  major  ions  NO+ 

and  0+.  These  diffusion  equations  are  derived  by  substituting  the  momentum  equations  (4 
and  5)  into  the  continuity  equations  (2  and  3)  tatog  into  account  charge  neutrality  (6)  and 
charge  conservation  (7).  Next,  all  of  the  coupling  and  nonlinear  terms  in  the  two  diffusion 
equations  are  expanded  in  a  Taylor  time-series,  but  the  algebra  is  considerably  more 
involved  than  that  described  above  for  the  E-region  solution  because  of  the  density 
derivative  terms.  Appendix  A  gives  the  resulting  diffusion  equations  after  they  have  been 
"linearized  in  time,"  This  set  of  coupled,  second-order,  parabolic,  partial  differential 
equations  is  solved  with  a  completely  implicit  numerical  technique,  as  shown  symbolically 
in  Figure  8  [FREGION2].  This  numerical  solution  requires  taking  derivatives  [GRAD], 
setting  up  coefficients  [ADCOEF,  RECOEF,  CBCOEF],  and  applying  boundary 
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conditions  [TOPDEN]  before  the  double  tridiagonal  matrix  is  solved.  The  bottom 
boundary  conditions  for  the  F-region  solution  is  obtained  from  the  top  of  the  E-region. 

After  the  NO+  and  0+  density  profiles  are  calculated,  the  N2+  and  O2+  density  profiles  are 
calculated  assuming  they  are  minor  ions  in  chemical  equilibrium  (equation  1). 

Figure  9  shows  an  example  of  a  complete  density  solution  with  different  solution 
procedures  and  spatial  steps  in  the  E  and  F  regions.  The  spaces  in  the  profiles  are 
intentional  to  more  clearly  show  the  two  regions.  For  this  example,  the  altitude  separating 
the  E  and  F  regions  is  180  km. 


5.3.  Temperatures 

The  numerical  solution  for  Ti  is  straightforward,  since  a  local  thermal  balance  is 
assumed  at  all  altitudes  (equation  8).  In  this  case,  the  equation  for  Ti  is  algebraic. 

The  numerical  solution  for  Te  is  more  involved  because  the  electron  energy  equation  is 
a  nonlinear,  second  order,  partial  differential  equation  (i.e.,  equation  9).  To  improve 
numerical  stability,  the  altitude  domain  is  divided  into  a  lower  domain,  where  thermal 
conduction  is  negligible  and  a  local  thermal  balance  exists,  and  an  upper  domain,  where  the 
complete  energy  equation  is  solved.  In  the  lower  altitude  domain,  Te  is  obtained  from  the 
resulting  nonlinear  algebraic  equation  by  iterating  from  an  initial  guess.  In  the  upper 

altitude  domain,  equation  9  is  solved  by  introducing  the  variable  0  =  by  expanding 
all  nonlinear  terms  in  a  Taylor  time-series,  by  replacing  all  derivatives  with  finite  difference 
expressions,  and  by  inverting  the  resulting  matrix  with  a  standard  tridiagonal  technique  (cf. 
Schunk,  1988a).  TTie  bottom  boundary  condition  is  obtained  from  the  Te  value  at  the  top 
of  the  lower  altitude  domain.  At  the  top  boundary,  a  downward  thermal  conduction  flux  is 
specified. 

Again,  for  numerical  stability  reasons,  the  spatial  grid  on  which  the  electron  energy 
equation  is  solved  is  decoupled  from  the  spatial  grid  on  which  the  diffusion  equations  are 
solved.  This  is  for  two  basic  reasons.  First,  the  altitude  separating  the  E  and  F  regions, 
Zef,  may  not  coincide  with  the  altitude  up  to  which  a  local  thermal  balance  applies,  Zc. 
Also,  because  the  density  and  temperature  gradients  have  different  spatial  scales, 
computational  efficiency  dictates  different  step  sizes  for  the  diffusion  and  energy  equations. 

Figure  10  shows  a  typical  spatial  grid  for  the  density  code  and  four  possible  spatial 
grids  for  the  Te  code.  For  the  density  code,  a  20  km  step  is  used  in  the  F-region  and  a  4 
km  step  in  the  E-region.  For  the  Te  code,  the  first  example  is  for  a  20  km  step  at  all 
altitudes  and  Zc  below  Zep;  the  second  example  is  for  AZ  =  20  km  and  Zc  above  Zef;  the 
third  and  fourth  examples  are  for  AZ  =  4  km  and  Zc  below  and  above  Zef,  respectively. 
Finally,  Figure  11  shows  Te  and  Ti  profiles  for  a  case  where  a  4  km  spatial  step  is 
necessary  in  the  Te  code  at  altitudes  above  Zc  because  of  a  rapid  increase  in  Te  with 
altitude.  The  spatial  steps  in  the  corresponding  density  code  are  4  km  in  the  E-region  and 
20  km  in  the  F-region. 


5.4.  H+  Formulation 

The  IFM  has  a  relatively  simple  prescription  for  calculating  H+  densities.  The 
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algorithm  developed  is  based  on  information  obtained  from  both  satellite  data  and 
numerous  polar  wind  studies  (cf.  Schunk,  1988b;  and  references  therein).  In  constructing 
Ae  model,  account  was  taken  of  the  fact  that  below  about  500-700  km,  H+  is  a  minor  ion 
in  chemical  equilibrium  at  all  latitudes  and  longitudes.  Since  the  main  source  and  loss  of 
H+  stems  from  the  accidentaUy  resonant  charge  exchange  reaction, 


0*  +  H«H*-»-0 

the  chemical  equilibrium  expression  for  H+  becomes 

N(H*)  =  1.13(T.  /Ti)‘''N(H)N(0*)/N(0) 


(10) 


(11) 


Therefore,  at  the  altitudes  where  this  expression  is  valid,  H+  can  be  obtained  using  the 
MSIS  atmospheric  model  values  for  N(H),  N(0),  and  Tn  in  combination  with  the  IFM 

values  for  N(0+)  and  Ti.  At  altitudes  above  the  domain  where  chemical  equilibrium 
prevads,  different  prescriptions  are  used  in  different  latitude  regions.  The  latitude  regions 
are  below  55o,  55o  to  65o,  and  above  65o  magnetic  latitude. 

At  low-mid  latitudes  (<55o),  the  H+  altitude  distribution  changes  quickly  from  chemical 
equilibrium  to  dMi^ive  equilibriuni  However,  as  it  turns  out,  when  H+  is  a  minor  ion  the 
chemical  equilibrium  and  diffusive  equilibrium  scale  heights  are  nearly  identical. 

Therefore,  with  a  negligible  error,  the  H+  density  can  be  calculated  with  the  chemical 

equilibrium  formula  (11)  up  to  an  altitude  where  N(H+)  =  N(0).  Above  this  transition 

altitude,  the  H+  density  decreases  exponentially  according  to  the  following  diffusive 
equilibrium  formula. 


N(H")  =  NB(H*)exp[(zB  -z)/Hj  (12) 

Hi=k(T.-HTi)/M(H")g  (13) 

where  Nb(H+)  is  the  bottom  boundary  value  obtained  from  the  top  of  the  chemical 
equilibrium  domain.  However,  in  order  to  obtain  a  smooth  transition  from  the  growing  to 
the  decreasing  exponential  solutions,  they  are  combined  over  a  small  altitude  region  near 
the  transition  altitude. 

At  high  latitudes  (>65®),  the  geomagnetic  field  lines  are  'open'  and  ,  as  a  consequence, 
H+  is  in  a  continual  state  of  outflow.  Both  model  studies  and  measurements  indicate  that 
the  H+  escape  flux  is  saturated  and  that  H+  remains  a  minor  ion  to  high  altitudes  (=3000 

km).  In  such  a  situation,  Ae  H+  density  first  increases  with  altitude  in  the  low-altitude 
chemical  equilibrium  domain,  reaches  a  peak,  and  then  decreases  with  altitude  with  a  scale 

height  that  is  equal  to  the  0+  scale  height.  An  examination  of  a  number  of  polar  wind 
studies  indicates  that  the  H+  peak  density  is  approximately  2%  of  the  CH-  density  at  the  peak 
altitude  (cf.  Schunk,  1988b).  Therefore,  the  simple  prescription  used  in  the  IFM  to 

calculate  H+  density  profiles  at  high  latitudes  is  as  follows.  Chemical  equilibrium  is 
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assumed  up  to  an  altitude  where  N(H+)  =  0.02N(O+).  At  high  altitudes,  the  H+  scale 
height  is  set  equal  to  the  0+  scale  height, 

H(H*)  =  H(0*)  =  k(T,  +T,)  /  M(OOg  (14) 

The  high  and  low  altitude  H+  density  segments  are  then  connected  with  a  parabolic  profile 
over  a  50  km  transition  region.  It  should  be  noted  that  the  resulting  H+  density  profiles  are 
very  similar  to  those  obtained  in  model  studies  of  the  polar  wind. 

The  latitudinal  region  55°  -  65o  corresponds  to  the  light  ion  trough,  wherein  the  H+ 
density  undergoes  a  transition  from  diffusive  equilibrium  at  55  o  to  continual  outflow  at 
650.  At  a  fixed  altitude,  a  satellite  traversing  this  region  from  low  to  high  latitudes  typically 
measures  a  constant  O  density,  but  an  H+  density  that  decreases  logarithmically  due  to  the 
increasing  importance  of  H+  outflow  (see  Figure  12).  With  allowance  for  this  additional 

feature,  the  prescription  used  to  calculate  H+  density  profiles  in  the  55®  -  65o  lati^de 
domain  is  very  similar  to  that  used  at  high  latitudes.  Specifically,  the  chemical  equilibrium 

formula  (1 1)  is  used  up  to  an  altitude  where  N(H+)  =  (factor)  x  N(C)+).  The  factor  varies 
logarithmically  with  magnetic  latitude  firom  1.0  at  55®  to  0.02  at  65o.  This  basically 
determines  the  H+  peak  density  and  height  At  high  altitudes,  the  H+  scale  height  is  set 

equal  to  the  O  scale  height  Finally,  the  high  and  low  altimde  H+  density  segments  are 
connected  with  a  parabolic  profile  over  a  50  km  transition  region.  This  latter  step 

determines  the  exact  values  for  the  H+  peak  density  and  height 

Figure  13  shows  typical  H+  and  0+  density  profiles  obtained  from  the  IFM  at  different 
latitudes.  The  geophysical  parameters  are  for  summer  daytime  conditions,  and  the  latitudes 
are  along  the  noon  magnetic  meridian  at  54o,  60®,  64®,  and  68°  (from  left  to  right).  The 
latitudes  correspond  to  one  below  the  light  ion  trough,  two  in  the  light  ion  trough,  and  one 

at  high  latitudes.  These  latitudes  cover  all  three  domains  in  the  H+  formulation  and,  hence, 
the  H+  density  profiles  in  Figure  13  are  representative  of  what  is  obtained  firom  the  IFM. 

5.5.  Equatorial  Model  Implementation 

As  already  described,  the  equatorial  model  developed  by  Anderson  [1973]  and  recently 
validated  by  Preble  etal.  [1994]  is  the  basis  for  the  IFM  equatorial  F-tegion  model.  Its 
implementation  into  the  IFM  required  the  following  st^s: 

1 .  The  Anderson  model  had  to  be  upgraded  from  a  single  longitude  to  a  global  model. 

2 .  The  Anderson  model  time  sequencing  had  to  be  modified  such 

th^t  global  ionospheric  "snapshots"  could  be  acquired  every  1  hour  and,  if 
requested,  a  forecast  as  short  as  1  hour  could  be  run. 

3 .  An  entirely  new  capability  of  "hot  starting"  the  Anderson  model  had  to  be 

developed. 

4 .  This  required  the  development  of  a  module  that  aeates  the  hot  start  database,  if 

required. 
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5 .  An  E-region  had  to  be  added  to  the  Anderson  F-region  O-  model. 

The  equatorial  model  is  interfaced  with  IFM  through  a  single  call  to  subroutine  IFM_ 
EQU.f.  This  subroutine  carries  out  various  operations  which  are  shown  in  block  diagram 
form  in  Figure  14,  The  first  of  which  is  to  interface  to  the  run  conditions  already  defined 
by  the  IFM.SETUP  file  and  augmented  by  the  EQU.SETUP,  At  this  point  in  the 
sequence,  the  equatorial  module  can  return  (in  the  event  of  no  equatorial  region  being 
needed);  compute  a  cold  start;  or  begin  with  a  hot  start.  The  hot  start  is  a  set  of  ascii 
datasets,  one  per  longitude,  that  contain  initialization  of  coefficients  for  the  Anderson 
model.  These  are  more  extensive  than  just  ionospheric  parameters,  ie,  density, 
temperature,  etc.  They  need  to  be  generated  by  running  the  Anderson  model  in  the  GET_ 
HOT_START  mode  (See  Figure  14).  Once  a  hot  start  is  available,  (note  if  a  cold  start  is 
required  the  subroutine  GET_HOT_START  is  run  first  to  obtain  the  corresponding  hot 
start  conditions),  the  subroutine  EQU_F_REGION  is  run  (see  Rgure  14).  This  generates 

a  database  of  field  line  0+  densities  at  the  required  Latitudes,  MLTs  and  UTs.  However, 
this  database  has  an  uneven  altitude  distribution,  hence,  a  further  stage  of  interpolation  is 
required.  EQU_ALT_INT  carries  out  this  interpolation  which  is  a  linear  interpolation,  on 

logic  (O^-)  and  linear  altitude  (see  Figure  14).  These  altitude  profiles  are  then  transferred  to 
the  final  IFM  output  database  as  0+  densities. 

The  final  stage  of  the  equatorial  model  is  to  add  an  E-region  to  the  above  F-region.  In 
Figure  14  this  is  shown  as  a  call  to  subroutine  EQU_E_REGION.  This  subroutine  uses 
the  IFM  mid  and  high  latitude  E  region  formulation  and  is  run  in  the  same  manner  as 
already  described  in  the  preceding  sections  on  the  E-region  in  the  mid-high  latitude 

formulation.  This  model  then  adds  the  molecular  ions  to  the  0+  in  the  final  IFM  output 
database. 

The  merging  of  the  equatorial  model  output  with  the  high-mid  latitude  model  output 
occurs  as  the  data  is  written  to  the  IFM  output  database.  This  is  done  within  the  EQU_ 
ALT_INT.f  and  IFM_E_REGION.f  subroutines.  There  is  no  requirement  for  a  final 
merging  program  to  match  the  various  model  outputs  into  the  final  IFM  output  database. 
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AZ  =  20  km 


AZ  =  4km 


Figure  6. 


The  IFM  is  based  on 


erent  altitude  grid  spacings  in  the  E  and  F  regions. 
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ALTITUDE  PARAMETERS 


Figure  7.  Block  dia^am  showing  the  sequence  of  calculations  done  at  15-minute 
intervals.  These  calculations  provide  the  required  inputs,  associated  parameters,  and 
plasma  temperatures  needed  by  the  density  solver. 
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Numerical  Solution 


Figure  8.  Block  diagram  showing  the  sequence  of  subroutine  calls  needed  to  obtain  a 
numerical  solution  of  the  ion  diffusion  equations. 
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Figure  9.  Altitude  profiles  of  the  molecular  and  oxygen  ion  densities  for  daytime 
steady  state  conditions.  The  density  profiles  were  obtained  with  different  solution 
procedures  in  the  E  and  F  regions. 


DENSITY  and  Te  GRIDS  DECOUPLED 

Density  Te  Te  Te  Te 

Code  Code  Code  Code  Code 


Zef:  E  region  -  F  region  interface 

Zc:  Heating  =  Cooling  below  this  altitude 


Figure  10.  Schematic  diagram  showing  the  decoupling  of  the  spatial  grids  used  in  the 
density  and  electron  temperature  solvers. 
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AZ  =  4km 


357/10:00 


SEC/IFM 


year  1987,  f10.7  70.0,  f10.7a  70.0,  ap  4.0,  kp  1.0 
Trajectory  SN  1,  Step  20.  NmF2  0.937E+04.  HmF2  248. 
MLT  16.85,  Latitude  71.2 


Figure  11.  Altitude  profiles  of  the  ion  and  electron  temperatures  at  a  location  and  time 
when  a  4  km  spatial  step  was  needed  because  of  the  steep  electron  temperature  gradient. 
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MID-HIGH  LATITUDE  TRANSITION 


UT  11:22  :25  :28  :31  :35 

D.LAT.  (S)  40°  50°  60°  70°  80° 


Figure  12.  Plasma  ion  measurements  from  the  ISIS  satellite. 


2 


Figure  13.  H+  and  O  density  profiles  along  the  noon  magnetic  ni^  tidian  at  latitudes  of  54o,  60o,  64o,  and  680  (from  left  to 


START 
IFM  EQU 


Figure  14.  Block  diagram  showing  the  logical  steps  needed  to  run  the  equatorial  part  of 
the  IFM. 


29 


6.  IFM  VALIDATION 


The  current  version  of  the  IFM  has  been  validated  at  three  different  levels.  The  first 
involves  validating  the  physics  and  chemistry  upon  which  the  numerical  model  is  based. 
The  next  is  verifying  that  the  streamlined  model  is  consistent  with  the  comprehensive 
ionospheric  model.  TTie  last  validation  concerns  comparisons  of  model  predictions  with  the 
kind  of  data  that  will  be  available  at  the  Space  Forecast  Center.  In  the  subsections  that 
follow,  our  validation  efforts  in  these  three  areas  will  be  briefly  discussed. 


6. 1.  Physics  and  Chemistry 

The  physics  and  chemistry  contained  in  the  comprehensive  ionospheric  model  have 
been  validated  over  a  14  year  period  and  the  results  have  been  published  in  a  series  of 
papers.  These  papers  are  Usted  in  Appendix  C  and  for  each  paper  highlights  are  given  in 
buUit  form.  The  appendix  only  covers  the  validation  efforts  concerning  mid  and  high 
latitude  data.  However,  the  Anderson  model  of  the  equatorial  ionosphere  has  also  been 
validated  over  the  years  and  it  is  clear  that  the  physics  and  chemistry  contained  in  that 
model  are  correct  (cf.  Preble  et  al,  1994;  and  references  therein). 

With  regard  to  the  mid  and  high  latitude  domains.  Appendix  C  shows  Aat  numerous 
satellites  and  radars  have  been  involved  in  various  aspects  of  our  validation.  Over  the 
years,  we  have  used  data  from  five  satellites  (AE-C,  DMSP  F2  &  F4,  NOAA-6,  DE-2), 
four  incoherent  scatter  radars  (Millstone  HiU,  Clhatanika,  EISCAT,  Areci^),  two  coherent 
radars,  and  about  50  ionosondes  spread  around  the  world.  Note  that  41  ionosondes  were 
used  in  one  study  alone  and  that  12  Soviet  ionosondes  were  used  in  another  study.  Note 
also  that  are  validation  efforts  covered  different  solar  cycle  (maximum  to  minimum), 
seasonal  (winter,  equinox,  sununer),  and  magnetic  activity  (high,  medium,  low) 
conditions.  Diurnal  variations  were  studied  as  well  as  longitudinal  effects.  Finally,  note 
that  the  validations  involved  the  model  predictions  for  NmF2,  hmF2,  the  electron  densities 
at  different  altitudes,  the  ion  composition,  and  the  electron  and  ion  temperatures. 

The  main  conclusion  to  be  drawn  from  the  various  validation  studies  that  we,  and 
others,  have  conducted  over  the  years  is  that  critical  inputs  are  needed  in  each  of  the  three 
latituchnal  domain.  If  these  inputs  are  known,  the  ionospheric  model  predictions  are  very 
reliable.  The  critical  inputs  are  the  plasma  convection  and  particle  precipitation  patterns  at 
high  latitudes,  the  thermospheric  wind  at  mid-latitudes,  and  the  electrodynamic  drift  at  low 
latitudes. 


6.2.  Streamlined  Versus  Comprehensive  Model 

As  part  of  the  PRISM  development  effort,  our  comprehensive  ionospheric  model  was 
run  54  times  in  order  to  produce  background  ionospheres  in  the  event  that  real-time  data 
were  not  available  to  drive  PRISM.  The  54  cases  corresponded  to  3  levels  of  solar  activity, 
3  seasons,  3  levels  of  geomagnetic  activity,  and  2  background  plasma  convection  patterns 
(IMF  By  positive  and  negative).  The  54  simulations  displayed  a  wide  range  of  ionospheric 
features,  including  tongues  of  ionization,  polar  holes,  nightside  and  dayside  ionization 
troughs,  auroral  ionization  enhancements,  elevated  sunlit  regions,  and  ion  and  electron 
temperature  hot  spots.  The  simulations  also  showed  how  these  features  varied  with 
altitude,  latitude,  longitude  and  universal  time. 

In  verifying  that  the  streamlined  and  comprehensive  models  produce  the  same  results. 
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we  repeated  many  of  these  54  cases  and,  indeed,  obtained  very  similar  results.  An  exact 
agreement  was  not  expected  because  the  streamlined  IFM  model  was  run  with  the  MSIS 
wind  (rather  than  the  simple  prescription  used  in  the  original  54  simulations)  and  because 

the  solar  EUV  fluxes  and  0+-0  collision  frequency  were  modified  according  to 
recommendations  from  the  PRIMO  community.  Fi^es  15  to  18  show  four  representative 
'  IFM  simulations  for  the  following  geophysical  conditions: 


Case  1 


Case  2 


Case  3 


Case  4 


Solar  Minimum 
Winter 
Low  Activity 


Solar  Maximum 
Winter 
Low  Activity 


Solar  Maximum 
Winter 
High  Activity 

Solar  Minimum 
Winter 
High  Activity 


(Fio.7  =  70) 
(Day  357) 
(Kp  =  1) 


(Fio.7  =  210) 
(Day  357) 
(Kp=l) 


(Fio.7  =  210) 
(Day  357) 
(Kp  =  6) 


(Fio.7  =  70) 
(Day  357) 
(Kp  =  6) 


All  four  cases  are  for  the  Northern  Hemisphere  and  magnetic  latitudes  from  40o  to  the  pole. 
For  cases  1-3,  snapshots  of  N(0)  at  300  km,  N(NO+)  at  300  km,  and  hniF2  are  shown  at 
0  UT.  For  case  4,  NmF2,  hmF2,  Ne  at  140  km,  Ne  at  800  km,  Ti  at  500  km,  and  Te  at 
230  km,  are  shown  at  both  0700  and  1900  UT. 

For  solar  minimum,  winter,  and  low  magnetic  activity  (Figure  15a-c),  the 
simulation  displays  all  of  the  features  that  are  expected,  based  on  the  corresponding 

simulation  conducted  with  the  comprehensive  ionospheric  model.  At  300  km,  the  0+ 
density  distribution  exhibits  elevated  densities  on  the  dayside,  a  polar  hole,  elevated 

densities  in  the  nocturnal  auroral  oval,  and  a  morning  mid-latitude  trough.  The  NO+ 
density  distribution  displays  similar  (not  identical)  features,  but  at  300  km  the  values  are 
considerably  lower  than  those  of  0+,  as  expected  for  low  magnetic  activity.  At  mid¬ 
latitudes,  hinF2  is  depressed  on  the  dayside  and  elevated  on  the  nightside,  which  is  in 
agreement  with  the  expected  behavior  due  to  the  transpolar  thermospheric  wind.  If  the 
solar  cycle  level  is  changed  from  minimum  to  maximum  keeping  all  of  the  other  parameters 
the  same  (Figures  16a-c),  the  O  and  NO+  density  features  are  similar  but  the  density  levels 
are  higher,  which  is  again  what  is  expected.  Likewise,  the  hmF2  values  at  solar  maximum 
are  considerably  higher  than  those  at  solar  minimum. 

The  simulations  that  are  shown  in  Figures  17  and  18  were  run  with  a  Volland  (1978) 
asymmetric  convection  pattern  with  enhanced  flow  in  the  dusk  section.  Because  of  this 
asymmetry  and  because  of  the  large  convection  speeds  associated  with  Kp  =  6,  the 
resulting  ionospheric  features  are  significantly  different  than  those  obtained  for  Kp  =  1. 
Comparing  the  corresponding  solar  maximum  winter  cases  for  Kp  =  1  (Figures  16a-c)  and 
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Kp  =  6  (Figures  17a-c)  several  important  differences  are  evident.  First,  an  0+  hole  and  a 

corresponding  NCH  density  buildup  appear  in  the  dusk  sector  coincident  with  the  strong 
convection  cell.  Also  note  that  the  polar  region  extends  to  lower  latitudes  for  Kp  =  6, 
which  results  in  lower  dayside  densities.  Finally,  note  diat  both  high  (500  km)  and  low 
(150  km)  hmF2  values  are  coincident  with  the  strong  convection  cell  in  the  dusk  sector. 

The  Kp  =  6  results  for  solar  minimum  are  shown  in  Figures  18a-f,  where  snapshots  of 
NinF2,  hmF2,  Ne  at  140  km,  Ne  at  800  km,  Tiiand  Te  are  displayed  at  0700  and  1900  UT. 
These  times  correspond  to  the  extreme  locations  of  the  terminator  in  the  magnetic  reference 
frame.  The  basic  features  shown  in  these  figures  are  very  similar  to  those  obtained  in  a 
corresponding  run  of  our  comprehensive  ionospheric  model  and  the  underlying  physics  is 
straightforward.  The  strong  convection  speeds  in  the  dusk  convection  cell  result  in 
elevated  ion  temperatures  due  to  ion-neutral  fnctional  interactions.  The  enhanced  Ti’s,  in 

turn,  cause  a  depletion  in  0+  and  an  enhancement  in  NO+  because  of  the  temperature 
dependence  of  the  C)++N2  iE  NO++N  reaction  rate.  The  elevated  Ti's  also  produce 
enhanced  Te's  due  to  collisional  coupling  of  the  ions  and  electrons.  However,  the  various 
density  and  temperature  features  are  modulated  by  the  motion  of  the  terminator  in  the 
magnetic  frame. 


6.3.  Observational  Data 

Ionospheric  models  are  being  compared  with  observations  on  both  national  and 
international  levels  in  order  to  validate  the  models  in  application  environments.  As  a  result, 
well  defined  and  quality  controlled  observation  datasets  are  available  against  which  the  IFM 
can  be  compared.  One  of  these  comparison  efforts  is  being  carried  out  jointly  by  the  NSF- 
CEDAR-PRIMO  working  group  and  the  STEP- VIM  worlMg  group. 

PRIMP:  Problems  Related  to  Ionospheric  Modelling  and  Observations  is  a  community 
wide  Ionospheric  model  and  observation  comparison  initiative  within  the  NSF-CEDAR 
program.  'Hiese  comparisons  have  been  done  for  solar  maximum  and  minimum,  summer 
and  winter  for  quiet  geomagnetic  conditions.  The  observations  have  come  from  the 
American  sector  (east  coast)  chain  of  digisonde  stations.  Comparisons  of  the  IFM  model 
results  with  observations  from  Millstone  Hill  and  with  the  USU  TDIM  model  are 
presented.  The  IFM  agrees  very  well  with  both  the  observations  and  TDIM.  At  this  time 
most  of  the  difference  between  the  IFM  and  TDIM  results  lie  in  using  an  oversimplified 
wind,  which  produces  too  large  a  midlatitude  nighttime  maintenance.  Tbis  is  especially  the 
case  during  solar  maximum.  With  the  Hedin  90  wind  model,  this  will  be  corrected. 
Figures  19a  through  19d  show  the  Ninf2  comparisons  of  IFM  (solid  line),  TDIM  (long 
dashed  line),  and  Millstone  HiU  observations  (dotted  line). 

Summary  of  PRIMO  Validation  as  of  1  August  1994. 

Solar  Minimum  -  Summer  (Figure  19a) 

Both  NmF2  and  hmF2  are  in  very  good  agreement  with  the  observations.  Indeed,  the 

agreement  appears  to  be  better  than  with  the  TDIM. 

Solar  Minimum  -  >^nter  (Figi^e  19b) 

NmF2  through  daytime  and  early  evening  is  good,  but  the  predawn  secondary 
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maximum  is  not  well  reproduced.  HmF2  may  be  a  little  high  (20  km).  Comparable  to 
TDIM  results. 

Solar  Maximum  -  Summer  (Figure  19c) 

Daytime  NmF2  and  hmF2  are  too  high.  The  default  wind  has  no  downward  drift  on  the 
dayside.  Such  a  wind  would  lower  the  F-layer  and  reduce  NmF2.  With  a  lower  sunset 
NmF2,  an  automatically  improved  nighttime  NmF2  would  be  achieved.  The  Hedin-90 
wind  model,  which  will  be  delivered  with  the  IFM,  has  a  daytime  wind  and  a  lower 
nighttime  wind. 

Solar  Maximum  -  Winter  (Figure  19d) 

Both  daytime  and  nighttime  NmF2  and  hinF2  are  too  large.  Nighttime  wind  is  too 
strong.  Replacing  wind  with  Hedin-90  would  reduce  nighttime  wind  and  introduce  a 
daytime  wind  that  would  lower  the  layer  and  decrease  the  daytime  NmF2- 

VTM:  Validation  of  Ionospheric  Models,  is  an  international  initiative  designed  to 
validate  both  physical  and  empirical  ionospheric  models.  The  VIM  datasets  were  used  for 
the  solar  maximum  Millstone  Hill  data  shown  in  Figme  19c  and  d.  VIM  also  provides  data 
at  other  locations  and  is  representative  of  the  data  available  at  the  AF  Space  Forecast  Center 
to  be  used  in  driving  ionospheric  models. 

F.OTTATORIAL:  A  very  extensive  validation  of  the  Anderson  model  was  carried  out  by 
Preble  et  al,  [1994].  They  used  an  extensive,  high  resolution  dataset  from  the  incoherent 
scatter  radar  at  Jicamarca,  Peru.  The  first  version  of  the  IFM  is  able  to  reproduce  these 
results,  verifying  that  the  restructuring  of  the  Anderson  equatorial  model  within  the  IFM 
has  been  successful.  Figure  20  compares  the  IFM  (top  panel)  with  the  Preble  et  al.,  [1994] 

(bottom  panel)  0+  density.  These  show  the  diurnal  altitude  variation  of  0+  over  Jicamarca 
for  solar  medium  conditions  Fio.7  =  128.7.  Preble  et  al.,  [1994]  demonstrated  that  their 

model  0+  variations  agree  very  well  with  the  observed  densities,  hence  validating  that  the 
IFM  will  do  the  same. 
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Figure  16a.  Snapshot  of  the  0+  density  distribution  at  300  Ion  and  0  UT.  The 
simulation  represents  solar  maximum,  winter,  and  low  magnetic  activity  conditions. 
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Midnight 


Figure  16c.  Snapshot  of  the  hmF2  at  0  UT.  The  conditions  are  the  same  as  for  Figure 
16a. 


39 


Noon 


(y  Density  (cm'®) 


Altitude 

Year 

Day 


LfT 

f10.7 

f10.7a 


Ap 

Kp 


300 

1982 

357 

0 

210 

210 

200 

6 


Figure  17a.  Snapshot  of  the  0+  density  distribution  at  300  and  0  UT.  The 
conditions  represent  solar  maximum,  winter,  and  high  magnetic  activity. 
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Figure  18a.  Snapshots  of  the  NmF2  distribution  at  0700  UT  (top  panel)  and  1900  UT 
(bottom  panel).  The  conditions  represent  solar  minimum,  winter,  and  high  magnetic 
activity. 


43 


Noon 


Midnight 


H„F,  (km) 

Altitude 

0 

Year 

1987 

Day 

357 

UT 

7 

f10.7 

70 

f10.7a 

70 

Ap 

200 

Kp 

6 

Module 

Version 

FT! 

2 

FTE 

2* 

Wind 

100 

B 

1 

E 

1 

Oval 

1 

COEF 

PRIMO 

Noon 


Midnight 


H„F,  (km) 

Altitude 

0 

Year 

1987 

Day 

357 

UT 

19 

f10.7 

70 

f10.7a 

70 

Ap 

200 

Kp 

6 

Module 

Version 

FTI 

2 

FTE 

2* 

Wind 

100 

B 

1 

E 

1 

Oval 

1 

COEF 

PRIMO 

Figure  18b.  Snapshots  of  the  hmF2  distribution  at  0700  UT  (top  panel)  and  1900  UT 
(bottom  panel).  The  conditions  are  the  same  as  for  Figure  18a. 
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Figure  18c.  Snapshots  of  the  electron  density  distribution  at  140  km  at  0700  UT  (top 
panel)  and  1900  UT  (bottom  panel).  The  conditions  are  the  same  as  for  Figure  18a. 
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Figure  18d.  Snapshots  of  the  electron  density  distribution  at  800  km  at  0700  UT  (top 
panel)  and  1900  UT  (bottom  panel).  The  conditions  are  the  same  as  for  Figure  18a. 
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Figure  18e.  Snapshots  of  the  ion  temperature  distribution  at  500  km  at  0700  UT  (top 
panel)  and  1900  UT  (bottom  panel).  The  conditions  are  the  same  as  for  Figure  18a. 
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Figure  18f.  Snapshots  of  the  electron  temperature  distribution  at  232  km  at  0700  UT 
(top  panel)  and  1900  UT  (bottom  panel).  The  conditions  are  the  same  as  for  Figure  18a. 
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£]  USU,  tdim;  data/ttdm.00005 
Observed  Data:  data/tobs.OOOl : 
_o  IFM  data/tiii.00006 


Fri  Sep  1 0  1 1 :23:1 8  1 993  PRIMO  project 


Figure  19a.  Diurnal  variations  of  NinF2  for  Millstone  Hill  as  obtained  from  the  IFM 
(solid  line),  the  TDEVi  (long  dashed  line)  and  incoherent  scatter  radar  measurements  (dotted 
line).  The  conditions  are  for  solar  minimum  and  summer. 
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Figure  19b.  Same  as  Figure  19a  but  for  solar  minimum  and  winter  conditions. 
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Figure  19c.  Same  as  Figure  19a  but  for  solar  maximum  and  summer  conditions. 
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Figure  19d.  Same  as  Figure  19a  but  for  solar  maximum  and  winter  conditions. 
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Figure  20.  Comparison  of  the  diurnal  variation  of  the  O  density  over  Jicamarca  as 
calculated  by  the  IFM  (top  panel)  and  by  Preble  et  al  (1994)  using  the  Anderson  equatorial 
model  (bottom  panel).  The  results  are  for  medium  solar  conditions. 
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7.  OPERATING  MODE 


The  current  version  of  the  BFM  is  self-contained  and  can  be  run  with  the  following  few 
parameters: 

Kp  (from  3  hours  prior  to  the  start  time  through 

Ae  duration  of  the  forecast) 

Fio.7  (solar  flux) 

Day  (decimal  number) 

Year  (gives  solar  cycle  phase) 

Start  Time  (UT  hours  and  minutes) 

Forecast  Duration  (hours  and  minutes) 

The  model  can  also  be  run  with  either  a  cold  or  a  hot  start.  With  a  cold  start,  no  initial 
global  distribution  of  densities  is  required.  However,  if  ionospheric  densities  from  the 
PRISM  model  are  available  or  if  the  output  from  a  previous  IFM  run  is  available,  the  IFM 
can  start  from  these  densities  and  then  forecast  the  subsequent  change  in  the  global  density 
distribution.  The  initial  global  distribution  of  densities  must  be  specified  on  the  ICED  grid. 

The  IFM  was  constructed  in  a  modular  form  so  that  when  real-time  data  are  routinely 
available  from  the  DMSP  satellites  and  the  DISS  and  TISS  ^ound-based  networks,  they 
can  be  ingested  to  provide  better  plasma  convection  and  particle  precipitation  patterns  as 
weU  as  better  real-time  density  updates.  The  next  phase  of  IFM  development  will  probably 
contain  this  improvement,  depending  on  data  availability. 

The  global  output  database  for  each  forecast  time  step  is  comprised  of  the  following: 

5  physical  parameters:  O,  molecular  ion,  H+,  Ti,  and  Te. 

35  altitude  steps:  at  different  altitudes  from  90  to  800  km. 

60  latitude  steps:  at  2o  intervals  poleward  of  50®  and  at  5o  intervals  equatorward  of  50o. 

24  longitude  (MLT/LT)  steps:  at  15o  intervals. 

This  corresponds  to  a  total  of  252,000  parameters  per  forecast.  Assuming  each  is  a  2  byte 
word,  this  is  a  0.5  Mbyte  dataset  If  the  IFM  is  running  with  a  forecast  every  15  minutes, 
then  the  IFM  generates  2  Mbytes/hour.  Any  one  of  these  datasets  can  be  used  to  do  a  iujl 
restart  of  the  B^.  In  24  hours  the  database  increases  substantially,  48  Mbytes. 

The  database  can  be  significantly  compressed  for  other  applications,  examples  of  which 
are  the  following  schemes: 

(1)  Peak  electron  density  maps  (1x1x60  x  24  x  4  x  2=1 1,520  bytes/hr). 

(2)  Ne  along  satellite*  orbit  (1x1x60  x  2x1. 5  x  2  =  360  bytes/orbit). 
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(3)  Northern  high  latitude  TEC  maps  (1x1x20  x  24  x  4  x  2  =  3,840  bytes/hr). 
*  assuming  an  ionospheric  satellite  orbit 
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APPENDIX  A. 


When  the  coupled  diffusion  equations  for 
O  and  NCH  are  linearized  in  time,  the  equations 
take  the  following  form: 
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APPENDIX  B 


SUBROUTINE  DESCRIPTIONS 
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ModvilsL 

railing  Name: 
PUTpQSgl 

Module: 
railing  Name: 
Purpose: 


Module: 
railing  Name: 
Purpose: 


Module: 
railing  Name: 
Purpose: 


Module: 
railing  Name: 
Purpose: 


Coefficients 


ADCOEF 


A  routine  that  calculates  the  coefficients  that  appear  in  the  coupled 
finite  difference  equations  for  O  and  NCH. 


Trajectory  Follow  -  advance  grid 


ANOTHERFLOW 


Prior  to  following  the  trajectories  for  an  IFM  time  interval  the  IFM 
grid  is  advanced  in  MLT  by  the  time  interval.  This  assures  that 
corotational  flux  tubes  stay  at  the  midpoints  of  their  cells  and 
reduces  the  number  of  flux  tubes  needed. 


Atmospheric  Densities  and  Temperatures. 

ATMOS  I 

A  routine  that  calculates  the  atmospheric  densities  and  temperatures 
as  a  function  of  altitude  at  a  specified  location  and  time.  A  fast 
routine  that  is  based  on  cubic  spline  fits  to  MSIS  outputs.  MSIS-90 
model. 


Auroral  Ionization 


AURORA! 


A  routine  that  calculates  the  auroral  ionization  rates  as  a  function  of 
altitude  at  a  specified  location  and  time.  A  fast  routine  which  has 
scaleable  altitude  profile  shapes.  The  scaling  factor  being  dependent 
on  the  auroral  energy  flux. 


Boundary  Conditions 

lED 

A  routine  that  provides  the  Iowct  boundary  densities  and  the  ion 
fluxes  at  the  top  boundary.  Chemical  equtiibrium  at  bottom 
boundary  and  zero  escape  fluxes  at  top  boundary. 
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Module: 
Calling  Name; 
P»rpc>S£l 

Modwis; 

Calling  Name: 
PvffpOSfil 


Module: 
Calling  Name: 

PTOO-SSI 


Module: 
Callinp  Name: 
Purpose: 


Module: 
Callinp  Name: 
Purpose: 


Coefficients 


CBCOEF 


Coefficients  associated  with  the  double  tridiagonal  matrix  inversion. 


Trajectory  Follow  -  uniformity  of  flux  tube  distribution 

cellcheck"! 


Every,  typically  15  minutes  of  model  time,  this  routine  checks  the 
distribution  of  flux  tubes  in  the  ionospheric  grid.  New  flux  tubes 
are  created  when  grids  are  empty  or  old  flux  tubes  dropped  when 
too  many  flux  tubes  are  in  a  single  grid  cell. 


Trajectory  Follow  -  define  cell  coordinates. 

|CELLXXX)RD^ 

Computes  the  latitude  (DLAT)  and  Magnetic  ^al  T^  (MLT)  of 
the  center  point  of  a  specific  grid  cell  on  a  variable  grid  or  a  variable 
grid. 


Chemical  Reactions 


CHEM  I 

A  routine  that  calculates  the  ion  and  electron  chemical  reactions  at  a 
specified  location  and  time.  A  fast  routine  which  accounts  for  only 
the  dominant  chemical  reactions. 


Diffusion  Coefficients 
IDIFF 

A  routine  that  calculates  the  ion  diffusion  coefficients  ^  a  function 
of  altitude  at  a  specified  location  and  time.  A  fast  routine  that  only 
includes  the  dominant  collision  processes  and  ignores  ion 
temperature  anisotropies. 
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Module: 
railing  Name: 

PvtrpQSSl 


Module: 
rallinp  Name: 
Purpose: 


Module: 
Tailing  Name: 
Purpose: 


Module: 
Tallinp  Name: 
Purpose: 


E-region  ion  and  electron  densities 
EREGION 

A  routine  that  calculates  the  ion  and  electron  density  profiles  in  the  E- 
region  at  a  specified  location  and  time.  The  NO+,  O2+,  N2+  and  On¬ 
ions  are  allowed  to  have  comparable  densities.  Chemical 
equilibrium  is  assumed. 


Trajectory  Follow  -  create  new  trajcetories. 


FILLCELLS  ~| 


During  the  trajectory  following  phase  when  an  empty  cell  is 
detected,  this  routine  creates  a  new  trajectory  at  the  center  of  that 
cell.  It  also  follows  that  trajectory  backwards  to  the  previous  time 
step  to  initialize  its  density,  etc.  The  procedure  is  as  follows: 

Rnds  empty  cell,  obtains  new  trajectory  coordinates 
(CELLCOORDS),  follows  trajectory  backwards  one  IFM  time 
interval  to  obtain  density  profile  (FLOWBACK),  then  reverses 
direction  and  covects  to  the  present  time  (REVERSE),  and  adds  it  to 
the  trajectory  data  base. 


Trajectory  Follow  -  trajectory  distribution  printout 


Ffi^LffiST  I 

Creates  record  of  the  distribution  of  flux  tubes  after  each  lE^  time 
period.  A  reduced  diagnostic  ouq)ut  con^atible  with  an  operational 
IFM.  This  diagnostic  would  only  be  called  if  needed. 


Trajectory  Follow  -  follow  all  trajectories. 

I  FLOW 

At  each  time  step  this  routine  must  generate  the  trajectories  for  all 

plasma  flux  tubes.  The  program  calls  a  heirarchy  of  routines 
including  the  modules  representing  the  magnetospheric  convection 
electric  field  (FLOWV3).  Corotation  is  added  in  the  FLOWVl 
routine.  Generic  interface  to  a  variety  of  convection  electric  field 
models.  Also  has  the  capability  to  handle  a  variable  time  increment 
when  following  a  trajectory. 
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Module: 
railing  Name! 
Purpose: 


Module: 
Calling  Name: 

CUII22SS1 


Module: 
Calling  Name: 
Purpose: 


Module: 
Calling  Name: 
Purpose: 


Trajectory  Follow  -  follow  a  trajectory  backwards. 
FLOWBACK  I 


Given  the  final  location  of  a  trajectory,  it  follows  the  trajectory 
backwards  in  time  to  obtain  its  initial  location  at  the  prior  IFM  time 
interval. 


Trajectory  Follow  -  convection  model 


FLOWV3 


Routine  to  specify  the  magnetospheric  electric  potential  in  the 
ionosphere.  The  Heppner  and  Maynard  [1987]  convection  model  is 
called  by  FLOWV3.  For  this  initial  level  the  input  is  kp  with  the 
assumption  that  the  A  pattern  is  appropriate.  At  this  level  no 
dependence  upon  the  interplaneta^  magnetic  field  is  used.  The 
variation  of  the  cross  polar  potential  widi  kp  follows  closely  the 
formula:  cross  tail  potential  =  (20+  14*  kp)  KV. 


F-region  ion  and  electron  densities 
FREGIONI 


A  routine  that  calculates  the  ion  and  electron  density  profiles  in  the  F- 
region  at  a  specified  location  and  time.  is  assumed  to  be  the  only 
major  ion.  The  NO  formulation  includes  diffusion,  but  it  is 
assumed  to  be  a  major  ion.  O2+  and  N2+  are  assumed  to  be  major 
ions  in  chemical  equilibrium 


F-region  ion  and  electron  densities 
FREGIOI^ 


A  routine  that  calculates  the  ion  and  electron  density  profiles  in  the  F- 

region  at  a  specified  location  and  time.  Both  O  and  NO  are 
assumed  to  be  major  ions  and  both  are  based  on  a  diffusion 
formulation.  O2+  and  N2+  are  assumed  to  be  minor  ions  in 
chemical  equilibrium. 
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Module: 

Callinp  Name: 

Purposs; 

Module: 
Tallin^  Name: 
Purpose: 


Module: 

railing  Name: 
Purpose: 

Module: 
Tailing  Name: 
Purpose: 


Module; 

Calling  Name: 
Purpose: 


Density  and  Temperature  Gradients 


GRADS  I 


A  routine  that  calculates  electron  and  ion  density  gradients  for  a 
multi-step  spatial  grid. 


Trajectory  Follow  -  initialize  the  convection  electric  field. 


INTTACrV 


Initialize  the  convection  electric  field  routines.  Define  the  kp  and 
other  proxy  indices.  IFM  uses  this  routine  to  interact  with  the 
AWS  geophysical  indices.  A  time  series  of  the  indices  is  read  in 
and  the  program  interpolates  as  needed  to  drive  IFM  with  a  time 
varying  skp  and  other  indices. 


Trajectory  Follow  -  grid  initialization 


initgridI 

Initialize  the  IFM  grid.  Genwic  global  routine  which  is  driven  by  a 
latitude  and  longitude  grid  specification. 


Trajectory  Follow  -  initialize  trajectories. 


INTITRAJ 


At  the  start  of  an  IFM  run  a  grid  is  specified  by  INTTGRID  and 
geophysical  conditions  by  INTTACrV.  These  data  are  then  used  by 
this  routine  to  initialize  trajectories  to  be  followed.  Places  flux  tubes 
at  the  center  of  grid  cells.  Can  be  run  with  an  uneven  grid  system. 


Altitude  Parameters 


ONE 


A  routine  that  calls  the  otiier  subroutines  which  provide  parameters 
as  a  function  of  altitude  (CHEM,  DIFF,  etc.)  at  a  specified  location 
and  time.  The  routine  also  calculates  derivatives  and  coefficients  for 
the  PDE's.  A  new  routine  that  allows  for  different  spatial  step  sizes 
in  the  E  and  F  regions. 
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Module: 
Calling  Name: 
Purpose: 

Module: 
Calling  Name: 
Purpose: 


Module: 
Calling  Name: 
Purpose: 


Module: 
Calling  Name: 
Purpose: 


Module: 
Calling  Name: 
Purpose: 


Auroral  Precipitation 

To  define  at  a  given  location  and  time  the  auroral  energy  flux  and  the 
auroral  characteristic  energy.  Improved  Hardy  auroral  oval  model. 


Riotoionization 

IphqtqI 

A  routine  that  calculates  the  solar  EUV  ionization  rates  as  a  function 
of  altitude  at  a  specified  location  and  time.  Ionization  rates  are 

calculated  only  for  zenith  angles  less  than  90o,  and  ionization 
firequencies  are  used  for  altitudes  above  300  km. 


Coefficients 


RECQEF 


A  routine  that  calculates  the  coefficients  for  a  double  tridiagonal 
matrix  inversion. 


Night  Sector  Ionization 


IRESONANT 


A  routine  that  calculates  the  night-time  ionization  rates  as  a  function 
of  altitude  at  a  q)ecific  location  and  time.  A  fixed  set  of  ionization 
profiles,  mainly  associated  with  resonantly  scattered  Lyman  Alpha. 


Trajectory  Follow  -  reverses  a  backwards  trajectory 


I  REVERSE 


Reverses  the  order  of  the  trajectory  steps  from  FLOWB  ACK  and 
adds  this  trajectory  to  the  trajectory  database. 
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Module:  Solar  Spectrum 

Calling  Name:  |SOLAR 

Purpose:  Defines  the  solar  spectrum  on  a  given  day.  A  39  wavelength 

spectrum  that  has  been  adapted  to  1 1  bins  and  used  in  the  TDIM. 
Very  fast,  depends  upon  FI  0.7. 


Module: 
Calling  Name: 

Pmpq^ 


Matrix  Inversion 
I  SOLVER 

A  routine  that  inverts  a  tridiagnol  matrix  for  one  partial  differential 
equation. 


Module:  Electron  Temperature 

Calling  Name:  ItbI 


Purpose:  A  routine  that  calculates  the  electron  temperature  as  a  function  of 

altitude  at  a  specified  location  and  time.  A  fast  routine  that  contains 
thermal  conduction,  but  includes  only  the  dominant  heating  and 
cooling  processes. 


Module:  Ion  Temperature 

Calling  Name:  HD 

Purpose:  A  routine  that  calculates  the  ion  temperature  as  a  function  of  altitude 

at  a  specified  location  and  time.  A  fast  routine  that  calculates  the  ion 
temperature  assuming  local  heating  and  cooling  processes  dominate. 


Module: 
^-ailing  Name: 


Top  Boundary  Conditions 
TOPDEN 


Purpose:  A  routine  that  calculates  the  top  boundary  conditions  for  O  and 

NO.  The  ionosphere  is  allowed  to  breathe,  but  no  ions  are  allowed 
to  escape  to  the  magnetosphere. 
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Module: 
Calling  Name; 


Neutral  Wind 
IWIND  I 


Purpose:  A  routine  which  provides  horizontal  wind  speeds  at 

location  and  time.  Altitude  variations  are  included, 
model. 


specified 
MSIS  wind 
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APPENDIX  C. 


Below  is  a  summary  of  the  papers  in  which  the  physics  and  chemistry  conned  in  the 
ionospheric  model  have  been  validated.  Highlights  of  the  studies  are  listed  in  buUit  form. 

HIGH  LATmJDES 

1 .  Sojka  et  al.  {JGR,  86,  2206,  198 1) 

•  AE-C  satellite  data  at  300  km 

•  High  latitude  regions  where  a  given  molecular  ion  dominates 

•  Model  predicted  correct  molecular  ion 

2.  Sojka  et  al.  (JGR,  87,  1711,  1982) 

•  DMSP  F2  and  F4  satellites  at  800  km 

•  Quiet  magnetic  conditions 

•  Diurnal  dependence  of  Ne  at  high  latitudes 

•  Model  predicted  correct  modulation  and  phase 

3.  Sojka  et  al.  (JGR,  88, 7783,  1983) 

•  Millstone  Hill  Radar 

•  Convection  pattern  obtained  from  radar  over  24-hour  period 

•  Corresponds  to  'average'  convection  pattern 

•  Model  Ne  at  500  km  in  general  agreement  with  measurement 

4.  Murdin  et  al.  (PSS,  32,  47,  1984) 

•  Chatanika  Radar 

•  Quiet  and  disturbed  days 

•  Inferred  vertical  ExS  &  wind  components 

•  Comparison  of  diurnal  variation  of  Ne 

•  Good  agreement  with  data  for  quiet  day 

•  Qualitative  agreement  with  data  for  disturbed  day 

5 .  Sojka  et  al.  (PSS,  33, 1375,  1985) 

.  DMSP  F4  satellite  at  800  km 
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•  Dayside  mid-latitude  trough  in  southern  &  northern  hemispheres 

•  Longitudinal  dependence,  depth,  and  asymmetry  between  North  and 
SouSi  troughs  were  in  agreement  with  measurements 

6.  Rasmussen  et  al.  {JGR,  91,  6986,  1986) 

•  Simultaneous  Chatanika  &  Millstone  Hill  data 

•  Coherent  scatter  radar 

•  NOAA-6  satellite 

•  Convection  pattern  from  radar  data 

•  Auroral  precipitation  pattern  from  satellite  data 

•  Generally  good  agreement  between  modelled  and  observed  Ne 
(25%) 

7.  Rasmussen  et  al.  (JGR,  93,  1922, 1988) 

•  Continuation  of  study  6  for  Te  &  Ti 

•  .  Very  good  agreement  between  measured  &  modelled  temperatures 

•  Topside  heat  flux  is  important 

8 .  Berkey  et  al.  (Polar  Res.,  48, 278, 1987) 

•  Halley  Bay  &  Siple  HF  Sounders 

•  Summer  &  Winter  Conditions 

•  NmFz  &  hinF2  comparisons 

•  Generally  good  agreement  between  modelled  and  measured  values 

•  Summer  hmF2’s  sensitive  to  wind 

9.  Sica  et  al.  (JATP,  50, 141, 1988) 

•  EISCAT  radar 

•  Large  downward  ion  drifts  [>400  m/s]  parallel  to  B  that  persist  for 
hours  deduced  from  data 

•  Model  shows  measurements  were  probably  not  interpreted  correctly 

10.  Rasmussen  et  al.  {JGR,  93,  9831,  1988) 

•  Arecibo  and  Chatanika  Radars 
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•  Compared  modelled  and  measured  densities  in  85-220  km  range 

•  Solar  maximum  &  minimum 

•  Generally  good  agreement  for  a  wide  range  of  conditions 

1 1 .  Sojka  et  al  {JGR,  95,  1 5275,  1990) 

•  IGY  lonosondes 

•  Mid-latitude  trough  exhibits  strong  magnetic  activity  and  longitude 
variations  (Whalen,  1987) 

•  Very  good  agreement  between  modelled  and  measured  variations 

12.  Sojka  et  al.  {Adv.  Space  Res.,  11,  39, 1991) 

•  Numerous  satellites 

•  Hoegy  &  Grebowsky  ( 1989)  mapped  morphology  of  ionospheric 
polar  hole  for  different  solar  cycle,  seasonal  UT,  and  Kp  conditions 

•  Model  correctly  predicted  the  observed  features 

13.  Sojka  et  al.  {JGR,  97,  1245,  1992) 

•  DE-2  Satellite 

•  E-field  &  particle  data  along  orbit 

•  Auroral  image  sequence 

•  Convection  &  precipitation  inputs  to  model  constrained 

•  Modelled  &  measured  Ne  coii5)ared 

•  Generally  good  agreement  for  quiet  day 

•  Discrepancies  during  substorm  related  to  lack  of  adequate 
convection  input 

MID-LATITUDES 

14.  Sojka  et  al.  (JATP,  50, 1027, 1988) 

•  Argentine  Islands  lonosonde 

•  (Juiet  magnetic  activity 

•  Monthly  diurnal  foF2  curves  from  solar  minimum  to  maximum 
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•  Excellent  agreement  between  modelled  &  measured  foF2 

•  Neutral  wind  is  critical  for  modelling 

15.  Wilkinson  et  al.  (Ann/.  Geophys.,  6,  31, 1988) 

•  10  Low  and  Mid-Latitude  lonosondes 

•  NmF2  diurnal  variation 

•  Mid-latitude  comparisons  depend  critically  on  wind 

•  Low-latitude  comparisons  depend  critically  on  ExB  drifts 

16.  Sica  et  al.  (JGR,  95,  8271, 1990) 

•  41  lonosondes 

•  Mid-latitudes 

•  Simultaneous  data 

•  Northern  and  Southern  hemispheres 

•  Diurnal  variation  of  NinF2  compared 

•  Winds  are  critical  for  modelling 

17.  Szuszczewicz  et  al.  (Adv.  Space  Res.  in  press,  1993) 

•  12  Soviet  lonosondes 

•  30-day-averaged  diurnal  NmF2 

•  Excellent  agreement  with  IRI 

•  Winds  are  critical  for  modelling 
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